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We examine a length scale that characterizes the spatial extent of heterogeneous dynamics in 
a glass-forming binary hard-sphere mixture up to the mode-coupling volume fraction cj> c . First, 
we characterize the system's dynamics. Then, we utilize a new method [Phys. Rev. Lett. 105, 
217801 (2010)] to extract and analyze the ensemble independent dynamic susceptibility x*{t) an d the 
dynamic correlation length £(t) for a range of times between the j3 and a relaxation times. We find 
that in this time range the dynamic correlation length follows a volume fraction independent curve 
~ ln(i). For longer times, £(t) departs from this curve and remains constant up to the largest 
time at which we can determine the length accurately. In addition to the previously established 
correlation r a ~ exp[£(r a )] between the a relaxation time, r a , and the dynamic correlation length 
at this time, £(r a ), we also find a similar correlation for the diffusion coefficient D ~ exp[£(r a ) ] 
with 8 ~ 0.6. We discuss the relevance of these findings for different theories of the glass transition. 

PACS numbers: 61.20.Lc,61.20. Ja,64.70.P- 



I. INTRODUCTION 

It is becoming increasingly apparent that growing 
length scales can be associated with the dramatic slow- 
ing down of the dynamics in glass-forming systems. One 
such length scale characterizes the spatial extent of the 
so-called dynamic heterogeneity. It has been found that 
upon approaching the glass transition the particles' mo- 
tion becomes increasingly heterogeneous and the parti- 
cles can be divided into "slow" and "fast" sub-sets [THS]. 
These sub-sets can be seen as distinct peaks in the prob- 
ability of the logarithm of single particle displacements 
P[log 10 (<5r); t] [6 -8 . Importantly, the slow and fast par- 
ticles are not uniformly distributed in space, but form 
clusters whose size increases as the dynamics slows. The 
average spatial extent of the clusters of slow particles can 
be defined as a dynamic correlation length. This dynamic 
correlation length and other closely related lengths have 
been studied in simulations [9lfl6]. experiments [T7lU9| . 
and discussed theoretically 

One convenient way to characterize the spatial extent 
of the clusters is to identify the slow particles and then 
determine their spatial correlations. In simulational in- 
vestigations this is typically done by analyzing the so- 
called four-point dynamic structure factors S<i(q;t). The 
four-point structure factor quantifies spatial correlations 
between the slow particles. The label "four-point" refers 
to the fact that S4,(q;t) is a correlation function of two 
two-point functions that are used to characterize parti- 
cles' dynamics and to define slow particles. Examples 
of these two point functions are the microscopic inter- 
mediate scattering functions and overlap functions, but 
other functions have also been used in the literature. 
Of particular interest are two quantities that can be ex- 
pressed in terms of S4 (<?;£): the dynamic susceptibility, 
Xi{t) = lim^o Sn(q\ t), which is a measure of the overall 
strength of the dynamic heterogeneity, and the dynamic 
correlation length, £(t), which characterizes its spatial 



extent. The relationship between X4(t) an d £(t) provides 
insight into the fractal dimension of the slow particles 
clusters. 

In spite of a relatively straightforward definition of 
Xi{t), its direct simulational evaluation suffers from a 
technical difficulty. The difficulty originates from the 
fact that in a typical simulational ensemble some global 
fluctuations are suppressed. For example, in most simu- 
lations of glass-forming liquids the number of particles 
and the volume is kept constant, thus the density of 
the system is constant and global density fluctuations do 
not contribute to the direct simulational calculation of 
Xi{t)- Berthier et al. [17j proposed that these suppressed 
fluctuations can be calculated utilizing a procedure de- 
rived by Lebowitz et al. |26j . This procedure results 
in a two-part expression for the dynamic susceptibility, 
Xi{t) = X4.(t)\x + %(t) where X4(*)|x is the susceptibility 
in an ensemble with x kept fixed (which can readily be 
obtained from simulations) and X(t) is a correction term. 

Berthier et al. [17] furthermore noted that while 
X4(0lx cannot be easily determined experimentally, the 
correction term X (i) can. This fact, together with the 
positive definite character of the former term, X4(^)|x > 
0, provides an experimental lower bound for Xi(t)- How- 
ever, the relative size of Xi{t) an d t ne correction term 
remained an open question. In addition, even if the ex- 
perimental lower bound was a good estimate of X4(*)i 
there was no reliable correlation between Xi(t) an d the 
dynamic correlation length £ (t) . 

The relative size of the two terms contributing to Xi(t) 
was investigated by Berthier et al. [22, 23j and by Bram- 
billa et al. [2Z]- The main conclusion was that as the 
dynamics slows, the correction term becomes an increas- 
ingly better approximation for the ensemble independent 
Xi{t)- The assessment and extension of this result is one 
of the subjects of the present paper. 

There were several earlier simulational investigations 
[TUHT5] of the correlation between the dynamic suscepti- 
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bility Xiif) an< i the correlation length £(t) and between 
the average dynamics (as characterized by, e.g. the a 
relaxation time, r a ) and but their results, by and 
large, disagreed [2S]. Recently, it has been realized that 
in order to get reliable results for X4(*) an d £(*) one nas 
to simulate systems considerably larger than was custom- 
ary [ini HHj • In an earlier short note [53] i we described an 
application of the method of Lebowitz et al. [25] to facil- 
itate the determination of £(t q ) using large scale, 80 000 
particles, simulations. We found that £(t q ) ~ ln(r Q ) 
over the full range of densities studied. This slower, log- 
arithmic growth of £(t q ) with r Q is more consistent with 
experimental findings than the power law growth found 
in many previous simulations. 

In the present paper we give details omitted in Ref. [2"§] 
due to length restrictions. In addition, we analyze the 
time dependence of both the dynamic susceptibility and 
correlation length, and investigate additional correlations 
between the length and the average dynamics. 

We describe the system and simulation method in 
Sec. |ll| and we briefly characterize the system's dynamics 
in Sec. |III| Since there is no accepted theory of the glass 
transition, examination of the system's dynamics can get 
bogged down with an extensive number of different fits 
and fit parameters. Throughout much of the paper we 
refer to two regimes: a mode-coupling like regime where 
power laws describe the data well, and a different dy- 
namic regime where the mode-coupling like power laws 
are not applicable. We discuss this characterization of 
the data and describe the relevant fits in Appendix |A"| 
After describing the system's dynamics, in Sec. |IV| wc 
investigate the dynamic susceptibility Xi(t) an d the dy- 
namic correlation length The technical details of 
the calculation of Xi(t) an d £(t) a re given in Appendix 
[B] In Sec. [V] we explore the connections between the dy- 
namic correlation length and the average dynamics. We 
finish with a discussion in Sec. IVII 



II. SIMULATION DETAILS 

We simulated a system introduced by Brambilla et al. 
[27j : a 50:50 binary hard-sphere mixture where the diam- 
eter d,2 of the larger sphere is 1.4 times larger than the 
diameter d\ of the smaller sphere. The size difference 
is chosen to inhibit crystallization. We studied systems 
with N = (Ni + N 2 ) = 80 000 particles at volume frac- 
tions <t> = n(Nidl+N 2 ti%)/6V equal to 0.4, 0.45, 0.5, 0.52, 
0.55, 0.56, 0.57, 0.58, and 0.59 and systems with 10 000 
particles at volume fractions <f) equal to 0.54, 0.575, 0.58, 
0.585, and 0.59. Additional simulations were performed 
at slightly different volume fractions and concentrations 
to obtain the derivatives needed in this work. To de- 
termine the derivatives with respect to <j>, we performed 
simulations at <j>±8<j) where 8<f> = 0.001 for <fi < 0.58 and 
6<fi = 0.0005 for 4> > 0.585. To determine the derivatives 
with respect to concentration c = N\/N, we performed 
simulations at c = 0.5 ± 0.05 for 6 < 0.58. The con- 



centration derivatives had very little 4> dependence over 
the range we examined. We found that they were not 
necessary to obtain accurate correlation lengths and sus- 
ceptibilities for (j> > 0.56. Thus we did not determine the 
concentration derivatives for <f> > 0.585. 

We performed Monte Carlo simulations with the lo- 
cal trial displacements of particles randomly chosen from 
a cube of length 0.1<ii. It has been shown that Monte 
Carlo dynamics reproduces well the long time dynamics 
of glass forming systems [3D]. Moreover, Brambilla et al. 
[2 7) have shown that the present system with this par- 
ticular Monte Carlo dynamics reproduces well the long- 
time dynamics of their experimental system - a dense 
poly-disperse hard sphere system in which hydrodynamic 
interactions can be neglected. 

The simulations were run for at least 100t q (r a is de- 



fined in Section III ) after the systems stopped aging. To 



check for the presence of aging, we examined two point 
and four point quantities to see if they significantly de- 
pended on the initial time of the calculation. We found 
that the dynamic susceptibility Xi(t)\(/>,c is very sensitive 
to aging, thus providing a good test of equilibration. We 
ran at least four production runs at each volume fraction, 
and the results are an average over those runs. Results 
are presented in reduced units where the unit of length 
is di and the unit of time t is one Monte Carlo step (a 
Monte Carlo step is one attempted move per particle). 
Since the center of mass of the system can drift, all po- 
sitions are calculated with respect to the center of mass 
(for the calculation of the center of mass position masses 
of all the particle were taken as identical). 



III. SINGLE PARTICLE DYNAMICS 

In this section we examine the slowing down of the av- 
erage dynamics. In addition, we show that there is an 
indication of dynamic heterogeneity in two-point func- 
tions, and the dynamic heterogeneity is increasing with 
volume fraction. 

We start by examining the volume fraction dependence 
of the a relaxation time, r Q , determined by a charac- 
teristic decay time of an average overlap function. The 
average overlap function is defined as 




Fo(t) = 



where w n (t) is a microscopic overlap function, 
w n (t) = e[a-|r n (f)-r n (0)|]. 



(1) 



(2) 



Here 0(x) is Heaviside's step function and r n (t) is the 
position of particle n at a time t. The microscopic overlap 
function uu n (t) select particles that did not move farther 
than a from their original positions during the time t. In 
this work we use a = 0.3. Correspondingly, the average 
overlap function F a (t) measures the average fraction of 
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particles which did not move farther than a from their 
original positions during the time t. We will refer to 
particles which did not move farther than a during time 
t as the slow particles. Thus, 



JV 



(3) 



is the number of slow particles during time t, and 
(N s (t)) = NF Q (t) is the average number of slow parti- 
cles during time t. 

F (t) encodes similar information as the 
self intermediate scattering function F s (q;t) — 



N ~ X (Er 



-q-[r„(t)-r n (0)] 



). Thus, it displays simi- 



lar characteristics. At high densities a pleateau region 
develops in the time dependence of F Q (t). Moreover, an 
early /3 relaxation regime can be identified as the decay 
to the plateau, and then the late /3 regime can be seen 
as a decay from the plateau. The characteristic time of 
the final decay from the plateau is referred to as the a 
relaxation time r Q . We define r Q adopting the formula 
used before for the self-intermediate scattering function, 



F {r a ) 



Consequently, according to this definition 



the average fraction of slow particles during time r a is 
about 37%. 

Shown in Fig.Qa) is F a (t) for = 0.4, 0.45, 0.5, 0.52, 
0.54, 0.55, 0.56, 0.57, 0.575, 0.58, 0.585, and 0.59 listed 
from left to right. For small volume fractions the de- 
cay is nearly exponential. At higher volume fractions 
the long time decay follows a stretched exponential form 
cxp[— (t/r) 13 ] with a weakly (j> dependent j3 s» 0.55. The 
stretched exponential relaxation is usually interpreted as 
an indication of dynamic heterogeneity. 

The other common way to examine the average dy- 
namics is to investigate the mean square displacement 
displacement 





(4) 



FIG. 1: (a): The average overlap function F a (t) for <f> = 0.4, 
0.45, 0.5, 0.52, 0.54, 0.55, 0.56, 0.57, 0.575, 0.58, 0.585, and 
0.59 listed from left to right. The dashed lines are stretched 
exponential, exp[— (t/r) 13 }, fits at <j> = 0.58 (/3 = 0.56) and 
4> — 0.59 (/3 = 0.54). (b): The mean square displacement 
(5r 2 (t)) for 4> = 0.4, 0.45, 0.5, 0.52, 0.54, 0.55, 0.56, 0.57, 
0.575, 0.58, 0.585, and 0.59 listed from left to right. 



which is shown in Fig. TJb). Again, for large <f> a 
plateau forms at intermediate times, then at long times 
(5r 2 (i)) = 6Dt where D is the self diffusion coefficient. 
Both the previously mentioned plateau in the average 
overlap function F a (t) and the plateau in the mean square 
displacement (<5r 2 (i)) are associated with the so-called 
cage effect where particles are temporarily trapped by 
cages of neighboring particles. We use the long time 
limiting behavior of (<5r 2 (i)) to obtain the self diffusion 
coefficient D. We define the j3 relaxation time Tp as the 
inflection point of ln[(5r 2 (<))] versus hi(t) (we found that 
it is easier to determine Tp from the ln[(5r 2 (t))] inflection 
point rather than from the F a (t) inflection point). This 
inflection point could only be determined for (f> > 0.5. 

Shown in Fig. [2] is the volume fraction dependence 
of the relaxation time, r a , and the inverse of the self- 
diffusion coefficient, 1/D. As in many glass forming sys- 
tems, there is an range of volume fractions in which power 



laws provide good fits to the simulation data. Since the 
mode-coupling theory predicts power law divergences of 
both r a and 1/D [3T], the volume fraction where the 
power law fits diverge is referred to as the mode-coupling 
volume fraction <j> c . However, there is no true diver- 
gence at 0c, the mode-coupling transition is said to be 
avoided, and there is emergence of new behavior beyond 
4> c ■ To quantitatively identify a mode-coupling like re- 
gion of the dynamics we fit r a and 1/D to power laws 
a(4>c — 0)~ 7fT ' D > where 7 r and 7d denote the power law 
exponents for r Q and 1/D, respectively (see Appendix 
[A] for a detailed description of the fits). We found that 
power laws describe our data well for 0.55 < <fi < 0.58 
with (p c — 0.59. These fits are shown as dashed lines 
in Fig. [2j We also find that our results for the relax- 
ation time are consistent with a fit suggested by Berthier 
and Witten [35] and later used by Brambilla et al. |27j . 
t q = Too cxp[B/(4>o — 4>) 2 ]. This fit is shown as the solid 
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FIG. 2: The relaxation time r a and the inverse diffusion coef- 
ficient 1/D. The dashed lines are power law fits a((j> c — 0) -7 
and the solid line is a fit suggested by Berthier and Witten 
|32| . Too exp[B/((j>o — 4>) 2 ]- Inset: Stokes-Einstein violation: 
for small volume fractions D ~ whereas for higher vol- 
ume fraction D ~ t~ 0,65 . 



FIG. 3: The probability of the logarithm of single particle dis- 
placement for the small particles calculated at r a for several 
representative volume fractions. For larger volume fractions 
there appears a multi-peak structure of P[log 10 (Sr);t] which 
indicates the existence of sub-populations of slow and fast 
particles. 



line in Fig. [2] (again, see Appendix [A] for a detailed de- 
scription of this fit) 

In the inset in Fig. [2] we investigate the relation be- 
tween two quantities discussed above, the a relaxation 
time and the self-diffusion coefficient. For small </> we 
find that D ~ r^ 1 and thus the Stokes-Einstein rela- 
tion is obeyed. With increasing <fi there appears to be 
a crossover to a weaker dependence of the self-diffusion 
coefficient on the a relaxation time. Thus, the Stokes- 
Einstein relation is violated. The breakdown of this rela- 
tions is considered to be one of the hallmarks of dynamic 
heterogeneity. 

Quantitatively, we find that for large volume fractions 
D ~ T Q T cr where a « 0.65. We should note that for 
even larger (f> it may be found that a < 0.65, and our 
result should be considered an upper bound. A value 
of a = 0.77 was found in an experimental glass-former 
[55] . and kinetically-constrained lattice-gas models pre- 
dict a fragility dependent a with values between 0.58 and 
0.88 [34] , The Random-First-Order theory also predicts 
a fragility dependent a [55] , 

To further explore the heterogeneous dynamics we ex- 
amined the probability of the logarithm of single parti- 
cles displacements, P[\og 10 (Sr); t] at r Q . This probabil- 
ity distribution is related to the self van Hove correlation 
function, G s (Sr;t) = (Sr — [r„(0) — r„(i)]), through the 
relationship P[log 10 (Sr);t] = \n{10)4Tr5r 3 G s (Sr; t). The 
advantage of examining P[\og 10 (5r);t] is that for Fick- 
ian diffusion (i.e. for a Gaussian distribution of single 
particle displacements) its shape is independent of time 
and its peak value is constant and approximately equal 
to 2.13 [6 . Thus, the time-dependence of the shape of 
P[log 10 (5r); t] is clear evidence of non-Fickian motion. 
Furthermore, a multi-peak structure of P[log w (5r); t] is 



indication of the presence of distinct sub-populations of 
particles and, thus, heterogeneous dynamics. 

Shown in Fig. [3] is P[log 10 (Sr); r a ] calculated for the 
small particles at (j> = 0.5, 0.55, 0.57, 0.58, and 0.59. 
The behavior is similar to what has been observed in 
other simulated glass-formers [5HE1 13S] m that multiple 
peaks emerge close to and at <f> c . These peaks correspond 
to slow and fast particles. In Fig. [4] we show the time 
dependence of P[log 10 (Sr); t)] for the small particles (a) 
and large particles (b) for cf) = 0.59. Shown are times 
equal to O.Itq,, t q , 5t q , 10t q , and 50TQ,. The multiple 
peaks are evident for both types of particles. The peaks 
are less pronounced and occur at later times for the large 
particles. 

For long times we would expect the particles to un- 
dergo Fickian diffusion. To compare the measured 
P[log 10 (5r); r a ] with those corresponding to Fickian mo- 
tion we show in Fig. [3] the probability distributions cal- 
culated at 50r Q assuming Gaussian distributions of dis- 
placements with the same (5r 2 (i)). It is clear that while 
for the small particles the difference between the mea- 
sured distribution and the Fickian one is relatively small, 
a pronounced difference is observed for the large parti- 
cles. Thus, even at the relatively long time, 50T a , large 
particles' motion is significantly non-Fickian. It should 
be emphasized that this conclusion cannot be obtained 
by only investigating the time dependence of the mean 
square displacement which grows approximately linearly 
with time on this time scale. Finally, we recall that in an 
earlier study we showed that in a binary Lennard- Jones 
system the time scale associated with the onset of Fick- 
ian diffusion increases faster with decreasing temperature 
than the a relaxation time [37]. We expect that a cor- 
responding result, i.e. that the time scale for the onset 
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FIG. 4: The probability of the logarithm of single particle 
displacements for <f> = 0.59 at 0.1r a , T a , 5t q , 10r a , and 50r a 
listed from left to right. The thin solid line is P[log 10 (<5r); t] at 
50r a calculated for a Gaussian distribution of displacements 
with the same (<5r 2 ) as measured in the simulation, (a) small 
particles, (b) large particles. 



IV. DYNAMIC SUSCEPTIBILITY AND 
CORRELATION LENGTH 

To examine the spatial extent of the heterogeneous dy- 
namics, we start with a somewhat qualitative approach 
and look at clusters of slow particles during r a utilizing a 
somewhat arbitrary definition. To visualize these clusters 
we define the slow particles as those whose displacement 
|r(t) — r(0)| was less than a = 0.3 over a time t = r a . We 
then define two slow particles to be in the same cluster 
if their initial positions were less than d a p 4- A Q ^ apart 
where d a p = (d a +dp)/2 and we used A Qi g = 0.02. Shown 
in Fig. [5] are clusters of more than 20 slow particles for 
4> = 0.55 and <fi = 0.59. It is apparent that the slow 
particles form clusters; moreover, there are more large 
clusters at <f> = 0.59 than at <p = 0.55. 

The definition of the clusters shown in Fig. [5] is arbi- 
trary. Alternative definitions results in different clusters. 
For example, a more common definition uses the separa- 
tion of the initial positions of slow particles correspond- 
ing to the first minimum of the respective pair correlation 
function, 



9ap(r) = 




N a (Np 



(5) 



where V is the volume, r = |r|, r nm = r n — r m , and the 
sums are over particles of a and /? type. Using such a 
definition we find that the clusters span the entire sim- 
ulation box. This is not surprising since, by definition, 
during time r a on average 37% of the particles are slow. 
Thus within the first minimum of g a /3(f) of a given slow 
particle another slow particle is likely to be found. 

To examine clusters of slow particles somewhat more 
quantitatively one can generalize the pair correlation 
functions, Eq. ([5]), and define a correlation function in- 
volving slow particles only, 



G 4 (r;t) 



V 



(N s (t))((N s (t))-l) 



(6) 



of Fickian diffusion grows faster with increasing volume 
fraction than the a relaxation time, holds for the present 
system. 



While one sees clear indications of different sub- 
populations of slow and fast particles at the higher vol- 
ume fractions from the results presented in Figs. [3] and [4j 
one cannot determine how these slow and fast particles 
are distributed in space. In the next section we examine 
the spatial correlations amongst the slow particles. As we 
mentioned in the opening paragraph of the introduction, 
these particles form clusters and a dynamic correlation 
length can be associated with the average spatial extent 
of the clusters. 



w n (t)w m (t)5[r 



where microscopic single-particle overlap functions w n (t) 
select slow particles, and (N s (t)) is the average number of 
slow particles, see Eqs. ([2]|3|. Note that the summation 
111 Eq. (|| is over all, small and large, particles. 

The function G/±(r;t) is usually referred to as a four- 
point pair correlation function. Note that by definition, 
in the thermodynamic limit, G^(r;t) — > 1 as r — > 00. By 
examining G^(r;t) — 1 we can examine the correlations 
between slow particles. In particular, the spatial extent 
of these correlations manifests itself in a slower decay of 
Gi(r; t) — 1 for large r. 

Investigation of the extent of the slow particles corre- 
lations through a direct analysis of G4 (r) is complicated 
by finite size effects. In particular, in the finite system 
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are proportional to exp[r/^(r Q )]/r where are determined 
by fits to 5*4(9; T a)- 




FIG. 5: Slow particles' clusters at r Q containing more than 20 
particles identified using the algorithm described in the text 
at (f> = 0.55 (upper figure) and <f> — 0.59 (lower figure). The 
white spheres are the large particles, and the black spheres 
are the small particles. Particles not belonging to the clusters 
are shown as black dots. 



a quantitative examination of the slow particles correla- 
tions. For the latter task we found it more convenient to 
analyze the wave- vector dependent analog of G±(r]t). 

Shown in Fig. |] is G\{r\T a ) - 1 for = 0.57, 0.58, 
and 0.59. The slower decay for larger <fi is evident, which 
indicates a larger correlation length associated with the 
average size of the slow clusters with increasing <p. 

To quantitatively determine the spatial extent of corre- 
lations of the slow particles we examine the q dependence 
of the four-point structure factor, 



= N- 



N 



(W(q,t)W(-q,t))-\(W(q;t))f 



NV 



(7) 



where W(q;t) is the Fourier transform of the spatially 
resolved microscopic overlap function, 



W(q; t)=^2 w n(t) exp[-iq ■ r„(0)] 



(8) 



and Gi(q;t) is the Fourier transform of G±{r;t) — 1. 

In the following two sub-sections we discuss two quan- 
tities that can be obtained from the four-point structure 
factor: the dynamic susceptibility Xi(t) an d the dynamic 
correlation length £(t). 



canonical ensemble the limiting large r value of Gi(r;t) 
differs from 1 by a term inversely proportional to the sys- 
tem size. To correct for this effect in a somewhat quan- 
titative way we determine the large r limit of G±(r; t) by 
finding the average value of G±(r;t) — 1 from r = 25.5 
to half the box length and then subtract this average 
from G\(r;t). The four-point function corrected in this 
way is denoted by G%{r]t). This function is shown in 
Fig. [6] We should emphasize that unlike in some other 
studies [HBTHO] we do not use this four-point function for 



A. Dynamic Susceptibility Xi(t) 

The dynamic susceptibility, X4(^)i i s defined as the q 
limit of the four-point structure factor, 



Xi{t) = lim S±{q;t). 



(9) 



Since as q -> 0, W(q; t) -> J2 n w n(t) = N s (t) (note that 
in the preliminary report [H] we used W(t) to denote the 
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^ n UJ n (i)) the dynamic susceptibility measures the ratio 
of the fluctuation of the number of slow particles to the 
total number of particles, and we could formally write 



X4(t) = N-U(N*(t))-(N s (t)y 



(10) 



It should be emphasized that the right-hand-side of 
Eq. ( 10 1 depends on the ensemble. In the ensemble used 



in our study the number of particles of both species are 
kept constant or, alternatively, the volume fraction </> and 
the concentration c are kept constant. Thus, hereafter we 
will denote the right-hand-side of Eq. (10) by X4(*)U,o 



X^tUc^N- 1 ((Nt(t))-(N s (t)y 



(11) 



where it is implicitly understood that the ensemble at 
the right-hand-side is the constant </> and c simulational 
ensemble. It should be noted that while Xi{t)\4>,c can 
be easily calculated in a simulation, in order to deter- 
mine the ensemble-independent susceptibility Xi(i) one 
needs to perform a rather delicate extrapolation proce- 
dure lim g _,.o <S4ta t). 

We note here that the difference between X4(0U,c an d 
Xi{t) is t ne reciprocal space manifestation of the finite 
size and ensemble dependencies of the large r limit of 
the four-point correlation function Gi{r;i). 

Berthier et al. [T7J pointed out that the susceptibility 
Xi(t) can be determined without extrapolating Si(q;t) 
by using the formalism introduced in Ref. [26] . This 
procedure results in the following expression 

X4(t) = XA{t)U,c + xl(t)Hi+X^)Xc{t)H 2 + xl(t)H z 
+F 2 Q (t)H i + F (t) X< p(t)H 5 + F (t) Xc (t)H 6 (12) 

where xS) = dF {t)/d(j) and X c(t) = dF a (t)/dc. The 
volume fraction dependent, but time independent, quan- 
tities H n are linear functions of the partial structure fac- 
tors S a p{q) extrapolated to q = 0. Note that we changed 
notation from previous work, Ref. [29], where we used 
G n instead of 7f„.This was done to avoid confusion with 
the four-point correlation function Gi{r\t). 

In Appendix [B] we present a derivation of Eq. |12| give 
the explicit formulae for the quantities H n , and describe 
how these quantities were evaluated. In the same ap- 
pendix we also describe an extrapolation procedure that 
confirmed the consistency of the definition ^ and the 
expression ( 12 1. 

It was further argued by Berthier et al. |17j that 
Eq. ( 12 1 could be used to establish an experimental lower 

> the correction terms 
constitute a lower 



bound for x<t{t)- Since Xi{t)\<t>,c 
at the right-hand-side of Eq. \12\ 
bound for Xi(t)- Furthermore, since around the a re- 
laxation time the first correction term is the dominant 
one, we could neglect all the other terms and thus arrive 
at 




FIG. 7: Time dependence of the terms that contribute to 
Xi(t) for 4> = 0.57. The arrow indicate the /3 relaxation time 
T0 and the a relaxation time Ta- 



li it could be shown that the x\(p)Hx term dominates 
close to the glass transition, one would have a simple ap- 
proximation for the dynamic susceptibility. We start by 
examining the time and <j) dependence of the terms on the 



right hand side of Eq. ( 12 ) to examine this approximation 
in detail. 

Shown in Fi g.[7 jis Xi($)\<f>,c and all the correction terms 
given in Eq. ( |12| for a representative volume fraction 
4> = 0.57. For very short times, the F$ (t)i?4 term is the 
largest, as it must be since it is the only term not equal 
to zero at t = 0, but it monotonically decays to zero. By 
the (3 relaxation time, X4WU,c and the x^W-Hi term are 
the largest, and by the a relaxation time they are around 
an order of magnitude larger than the other terms. For 
the volume fraction shown, these two terms are almost 
equal around r Q , but the x^W-ffi term becomes larger 
at later times. 

We note that with increasing volume fraction 
xl(r a )H 1 grows faster than Xi(Ta)U,c- The X%{ T o)Hi 
term becomes the dominant one for (j> > 0.58, Fig. [8] 
This is qualitatively consistent with result of Brambilla 
et al. 2Z] • The quantitative difference between our Fig. [8] 
and results shown in Fig. 3a of Ref. [57] originates from 
the fact that Brambilla et al. systematically overesti- 
mated the isothermal compressibility which enters into 
their correction term. 

For our range of volume fractions we do not find that 
Xi( T a)\ct>,c is negligible compared to x^,( T a)F[i. However, 
for volume fractions larger than the ones examined in 
this study, it is likely that X4( r o) is we H approximated 
by X 2 ^a)Hi term alone. 

As can be seen in Fig. [jj around r a a good approxima- 
tion to X4W is 



xl{t) = xS)Wc + xl{t)H 1 . 



(14) 



Xi(t) > xl(t)Hi 



(13) 



The time dependence of this quantity is shown in Fig. [9] 
We observe that xfW grows with time, indicating an 





FIG. 8: Volume fraction dependence of the constant <j> and 
c part of the dynamic susceptibility, X4( t q)U,c (circles), the 
dominant correction term, X%{ T a)Hi (triangles), and the en- 
semble independent susceptibility x±{p) calculated from Eq. 
(12 l (open squares). 



FIG. 9: The approximation of the dynamic susceptibility 
XS(*) = X4(t)kc + X*(*)fli for $ = 0.5, 0.52, 0.54, 0.55, 0.56, 
0.57, 0.575, 0.58, 0.585, and 0.59 listed from left to right. This 
approximation is accurate around r a , i.e. around the peak 
shown in the figure, and becomes increasingly more accurate 
as 6 increases. 



increase of the overall strength of dynamic heterogeneity, 
until it reaches a peak that occurs around r a and then 
decreases to zero at later times. The decrease in xtW 
represents a diminishing of the overall strength of the 
heterogeneous dynamics but the length scale associated 
with slow clusters do not have to follow the same trends. 

Toninelli et al. [UJ and Chandler et al. [42] examined 
theoretical predictions for the time dependence of Xi(t) 
and compared them with, inter alia particle-based sim- 
ulations. It is not clear whether the latter comparisons 
were hindered by fact that global fluctuations were sup- 
pressed in simulations. However, in general, a common 
feature predicted by many theories is a power law growth 
of Xiif) while approaching the peak. This fact prompted 
us to look for power laws in xt(*)- 

For smaller volume fractions we do not find any region 
of power law growth approaching r a . Around <f> = 0.56 
there emerges a region where xfM appears to grow ac- 
cording to a power law, but the exponent in the power 
law depends on </>. This is due to the two contributions 
to x%(t) having different magnitudes and time dependen- 
cies. For example, for <j> = 0.59 we find that xl(*) ~ i ' 665 
in the a relaxation regime. This growth is due to a com- 
bination of x%(t) ~ *°' 75 an d Xi\<t>,c ~ t ' 55 m the a 
relaxation time regime. This analysis suggests that the 
power law growth of xi{t) does not necessarily have a 
deeper meaning, at least for volume fractions accessible 
in our study. 

As we remarked above, we expect that for sufficiently 
high volume fractions the growth of X±{t) can be ob- 
tained from experiments using the xj^W-Hi correction 
term as an approximation for Xi(t)- Moreover, if time- 
temperature superposition holds, then the growth of 
x|(t) is related to the growth of the a relaxation time. 
Below, we investigate the consequencies of this idea. 



For hard sphere systems, time-temperature superpo- 
sition is replaced by time-volume fraction superposition. 
Specifically, the statement is that F (t/r a ) overlaps in the 
a relaxation regime when plotted for different <f). We find 
good overlap for <f) > 0.58, although we do observe small 
systematic deviations. Moreover, we find that F (t/T a ) is 
well described by a stretched exponential Ae~ i ^l Tc '^ (see 
Fig. [I). Thus, ignoring the weak volume fraction depen- 
dence of A and f3 the first correction term in Eq. ( 12 1 at 
r Q is given by 



xl{r a )H x = 



2A/3 



dln(r a ) 



e- 2 #i. 



(15) 



Recall that this is the largest term for <j) > 0.58. Since 
liquid structure is weakly volume- fraction dependent, 
Hi changes slowly with <f>. Notice that A, (3 and Hi 
are all less than one or equal to one for all ef>, and 
are all less than one for <j> ^ 0.58. Note also that 
2/{d\ + df) = 0.53419/g^. Thus, the coefficient multi- 
plying [91n(r Q )/91n((/))] 2 is less than one at all </> and 
is very weakly 4> dependent. Consequently, Xi(t) be- 
haves as [91n(r Q )/91n(0)] 2 when the x% term is domi- 
nant. Finally, since at the largest volume fractions r Q = 
Too exp[j3/(</>o — 4>) 2 ] provides a good fir to our data, then 
these arguments indicate that Xi( T a) ~ 4 >2 ( ( t ) o ~ ( t>)^ & 
close to <pQ. We find that this indeed provides a good 
description of our results for the dynamic susceptibility 
(See Fig.fll]). 



B. Dynamic Correlation Length £(t) 

To define the dynamic correlation length £ (t) we need 
to examine the long wavelength (small wave-vector) be- 
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0.58. 



10: Dynamic structure factor Si(q;T a 
0.57, 0.56, 0.55, 0.55, 0.52, and 0.50. 



i for <f> = 0.59, 
The values at 



q = were calculated using Eq. ( 12 1. The lines are Ornstein- 
Zernicke fits. 



havior of the four-point structure factor. Specifically, we 
use the following definition of the dynamic correlation 
length: 



lim q 



_ 2 / linifc^o S 4 (k;t) 



1 



(16) 



Definition ( 16 ) is consistent with asymptotic small wave- 
vector Ornstein-Zernicke form of the four-point structure 
factor, 



Si(q;t) 



X4(t) 



as q 



0. 



(17) 



We should note that formally definition ( 16 1 is equivalent 
to defining £(i) as the second moment of Gi(r\t) — 1 
divided by the zeroth moment of G/±(r\t) — 1. However, 
we found that the finite size and ensemble effects are 
easier to account for in the reciprocal space and therefore 
we used definition (16). 



To obtain reliable results for X4(t) an d £(t) we fitted 
the simulation results to several different functional forms 
(see Appendix [B] for details) . We determined that the 
best procedure is to fit Si(q; t), including as q = value 
the right-hand-side of Eq. (12), to the Ornstein-Zernicke 



form while restricting the fitting range to q < 1.5/£(t). 
Such fits at r a are shown in Fig. [Toj for volume fractions 
4> = 0.59, 0.58, 0.57, 0.56, 0.55, 0.52, and 0.50, calculated 
using the 80 000 particle simulations (note the the system 
size dictates the smallest non-zero wave- vector) . 



Shown in Fig. 11 is the volume fraction dependence 
of the resulting dynamic susceptibility and correlation 
length. Note that in this figure we also included results 
obtained applying the same procedure to data obtained 
from 10 000 particle simulations. The consistency of both 
sets of results indicates that our procedure can be used 
to determine the dynamic susceptibility and correlation 
length using moderately large systems. 



We find that the volume fraction dependence of £(r a ) 
can be well described by many different fit functions (the 
results of the mode-coupling like fits are described in de- 
tail in Appendix |Aj. To be consistent with previous fits 
to t q and D, we fix <j) c — 0.59 and fit £(t q ) to a mode- 
coupling like power law, £(r Q ) ~ (4>c — 0)~ 7<: • This results 
in 7£ = 0.5 ± 0.1. The corresponding fit is shown as the 
dashed line in Fig. 11 The value of 75 obtained from the 
fit does not agree with the inhomogeneous mode-coupling 
prediction of 7^ = 0.25 [25 ] . 

Next we fit £(r Q ) to £0 + C(<po — 4>)~ 2 over the whole 
range of <p, which gives <j)o = 0.0635 ± 0.004 and £ = 
0.37 ± 0.1. The corresponding fit is shown as the solid 
line in Fig. 11 and provides an accurate description of 
£(T a ) for every volume fraction examined in this work. 
Note that the same <j>o = 0.635 was obtained from fits 
of the a relaxation time to the formula suggested by 
Bcrthier and Witten [32], r a = Too exp[i?(0o~ ( / , ) _2 )] (see 
Appendix [AJ. This observation suggests that the follow- 
ing correlation between the a relaxation time and the 
length, r a = To exp[fc£(r Q )]. We discuss this relationship 
in Sec. E] 

We now look at the scaling relationship between 
the dynamic susceptibility and the length, X4( T a) ~ 
?( T a) 2-,? - For compact clusters it is expected that 
2 — r\ = d where d is the spatial dimension. Shown in 
the inset to Fig. 11 is the scaling fit for <fi > 0.56. We 
obtain 2 — 7/ = 2.9 ± 0.1, which indeed suggests com- 
pact clusters. The comparison of this result with the 
inhomogeneous mode-coupling theory prediction |211 125] 
is a little involved. The theory analyzes a three-point 
susceptibility and finds that in the a relaxation regime 
limq^o % 3 (q; r Q ) ~ £(r a ) 4 . Field-theoretical considera- 
tions [S3] [33] indicate that the dynamic susceptibility is 
a quadratic function of the three-point susceptibility. A 
combination of both results would suggest a prediction 
Xi( T a) ~ £( T a) 8 which is clearly well outside the simula- 
tional result. 

Since we find that £(r Q ) = £ + C(0o ~ < f ) )~ 2 provides 
a good description of all the data and that Xi{ T a) ~ 
a4^(r Q ) 3 for 4> > 0.55, we show a 4 [£ + C((f> — 4>)~ 2 ] 3 
as the solid line through Xi{ T a) in Fig. [TTJ Note that 
this is not an independent fit. However, it describes the 
data fairly well over the whole range of studied volume 
fractions. Moreover, it is consistent with the limiting 
behavior Xiija) ~ <^ 2 (^n — <t>)~ 6 obtained from the first 



correction in Eq. (12) 



In Fig. [TT] we also show, as a dashed line, the third 
power of the mode coupling fit. As expected, it gives a 
reasonable description of the data between 0.55 < <fi < 
0.58. 

We now examine the time dependence of the dynamic 
correlation length. Shown in Fig. 12 is £(t) versus time 
for different (f>. For all (f), the dynamic correlation length 
grows with time and then plateaus at later times, and 
remains constant up to the maximum time at which we 
can evaluate it. We cannot accurately calculate £(i) for 
t > 10t q because there are few slow particles at such long 
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FIG. 11: The dynamic correlation length £(r a ) versus 6 (left 
axis) and the dynamic susceptibility Xi( T a) versus <j> (right 
axis). The solid line through £(r a ) is a fit to £(r a ) = £o + 
A(4>q — 4>)~ 2 where (j>o = 0.635. The dashed line through 
£(r a ) is a mode-coupling like fit ^(r a ) ~ (6 C — 0) _7? where 
4>c ~ 0.59 is fixed and we obtain 75 = 0.5±0.1. In the inset we 
show Xi( T a) versus £(r a ), and we find that Xi( T a) ~ £( r a) 2 ' 9 
for A > 0.55. 



times. Thus, we do not know the fate £ at later times. 

We note that one previous simulational investigation 
of the time dependence of the dynamic correlation length 
resulted in length whose time dependence roughly fol- 
lows that of the dynamic susceptibility [TO]. A different 
study found a monotonically increasing dynamic corre- 
lation length [41]. On the other hand, an earlier study 
[35], which used a somewhat different definition of the 
dynamic correlation length, found that the length was 
increasing with time but then plateaued after t q . A 
similar behavior was found in a very recent study of a 
two-dimensional lattice gas glassy system [ID] . The two 
latter results are (at least qualitatively) consistent with 
our findings. 

We remark that a plateau in the time dependence of 
a characteristic dynamic length is predicted by the in- 
homogeneous mode-coupling theory |21l 125] . However, 
the plateau predicted by this theory occurs around the 
/3 relaxation time, and not after the a relaxation time as 
seen here. 

There are two somewhat surprising features in the re- 
sults shown in Fig. [12] First, for (j> > 0.56 and between 
rg and r a , the dynamic correlation length is independent 
of the volume fraction and only depends on time. We 
don't have sufficient data for times smaller than the /? 
relaxation time, but we expect that this universal behav- 
ior breaks down at short times. The correlation length 
appears to follow a master curve until it reaches a vol- 
ume fraction dependent asymptotic value, which we will 
refer to £ max . We find that this master cure can be well 
described by £(i) = a\n(bt), and this fit is shown as a 
solid line in Fig. |l"2") 



FIG. 12: The dynamic correlation length versus time for & > 
0.52. The correlation length appears to follow a universal 
curve until it reaches a volume fraction dependent maximum 
value, and then it stays approximately constant at later times. 
The straight line is a fit to the data, £(i) = aln(bt). 



Second, we find that the time at which £(t) saturates, 
which we denote as r max , exceeds the a relaxation time 
and the time at which Xi{t) peaks, r pea k (we find that 
T a and Tpeak have the same volume fraction dependence) . 
Thus, the dynamic correlation length seems to be grow- 
ing further while the overall strength of dynamic het- 
erogeneity, measured by the susceptibility, is decreasing. 
This somewhat surprising finding means that, while at 
longer time scales there are few slow particles, the char- 
acteristic size of the clusters of these particles seems to 
be constant (at least up to 10r a ). 

Since r max exceeds r a , it is obvious that £ max exceeds 
£(tq,). Interestingly, there seems to be a linear rela- 
tionship between these two lengths, Fig. [l3j Combin- 
ing the linear relationship between £ max and £(tq,) with 
the previous observation that £,(t) = a\n(bt), we see that 
T"m ax ~ ''"a and through the fits of £ max we determine 
e = 1.3 ± 0.1. 

The time dependence of the dynamic correlation length 
shown in Fig. [12] suggests that possibly we should focus 
more on the plateau value of the dynamic correlation 
length, £ max , than on the length at the a relaxation time, 
£(t q ). This suggestion is left for future investigation. 



V. CORRELATION LENGTH AND AVERAGE 
DYNAMICS 

In this section we examine the relationships between 
the dynamic correlation length £,{r a ) and the two sim- 
plest quantities characterizing the average dynamics, the 
relaxation time, r a , and the self-diffusion coefficient, D. 
We note that most theoretical descriptions of glassy dy- 
namics focus on the temperature dependence of the dy- 
namics. Consequently, relationships between £(r Q ), r a 
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FIG. 13: The plateau value £ max versus £(t q ). The solid line 
is a linear fit. 



and D involve temperature. However, temperature is 
not a relevant control variable for our hard sphere sys- 
tem. Instead, the volume fraction is the usual control 
parameter and, therefore, in formulae discussed below 
we omitted temperature. 

We begin by examining the relationship between £(r a ) 
and t q . The mode-coupling theory predicts a power law 
T a ~ £( T a) z ■ In contrast, the Adam-Gibbs [33] and 
Random-First-Order- Transition theories [44] predict an 
exponential dependence of the relaxation time on a cor- 
relation length, £ s , describing the size of cooperatively 
rearranging regions, r a ~ exp(£f) where 9 = 3 in the 
Adam-Gibbs theory and is a parameter in the Random- 
First-Order- Transition theory. While it is currently un- 
clear if our dynamic correlation length is the same as 
the correlation length in Adam-Gibbs and Random-First- 
Order-Transition theories (in particular our length de- 
pends on time whereas £ s does not have an obvious time 
dependence) we examine relationships between £ and r a 
suggested by those theories. 

We find that a power law describes well the correlation 
between £ (r a ) and r Q over the mode-coupling regime with 
z = 4.8 ±0.3, see Fig. ~ 
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This exponent disagrees with 
the inhomogeneous mode-coupling theory prediction of 
z « 10 [2TJ [35] and with some of the previous simula- 
tional studies [101 02] , but it is consistent with other pre- 
vious simulational investigations QT| Q31 US] (note that 
in majority of earlier studies the inverses of the exponent 
z was given). 

We find, however, that an exponential dependence of 
r a on £(r a ) provides a better description of the data 
over a larger range of volume fractions. A fit to r Q = 
T exp(£; r £(T Q ) e ) gives 9 = 1.1 ±0.2. Thus, we fix 9 = 1.0 
and fit r Q = To exp(fc r £(r a )), which is shown as a solid 
line in Fig. 14 Note that the quality of the latter fit 
is fully consistent with the fact that independent fits 
r a = cxp[B(0 o - (j))- 2 ] and £(r a ) = £ + C((f> - </>)~ 2 
result in the same value 4>q — 0.635 (see Appendix [A| . 



FIG. 14: The alpha relaxation time r a and the diffusion co- 
efficient D as a function of £(r a ). The solid straight line is an 
exponential fits to the relaxation time data r a ~ exp[£(r a )] 
and the solid curved line is a fit 1/D ~ exp[£(ra) 9 ] with 
9 = 0.6. The dashed line is a mode-coupling theory like fit to 
t q ~ ^{T a ) z with z = 4.8. 



In Fig. 
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we also show 1/D versus £(T a ). The re- 
sults do not seem to follow a straight line and thus 
we do not find that D ~ exp(— fcrj£(T Q )), but rather 
D = D exp(-k D £(T a ) e ) where 9 = 0.61 ± 0.04. The 
latter fit is shown as a solid line in the figure. Again, we 
note that the quality of the self-diffusion coefficient fit 
is quite good. However, we shall also note that a com- 
bination of both correlations r a = tq exp(k T ^(r a )) and 
D = £>o exp(— fcu£(T Q ) ) is, strictly speaking, not com- 
patible with a power-law relationship between the self- 
diffusion coefficient and the relaxation time discussed in 
Sec. HI] 

Finally, we briefly mention two other, different inves- 
tigations that analyzed somewhat different characteristic 
dynamic lengths. 

Saltzman and Schweizer [5] HSJ 03] investigated a char- 
acteristic length associated with the onset of Fickian dif- 
fusion. They showed that this length, £d, depends log- 
arithmically on the relaxation time, £d ~ ln(r a ). This 
result is consistent with our relation between the dynamic 
correlation length at the a relaxation time and the a re- 
laxation time. 

A similar crossover length was examined in models of 
facilitated dynamics [47] . Berthier et al. defined a length 
scale associated with the onset of Fickian diffusion as 
I* ~ VDt^- They noted that t - T { a~ a)/2 where 
a is the dynamic exponent describing the violation of 
the Stokes-Einstein relation. Thus, for our system one 
would expect I* ~ t^ 2 . We do find that for our sys- 
tem the relation y/Dr a ~ t°' 2 is obeyed for large r a . 
However, the length £* has different volume fraction (or 
relaxation time) dependence from our dynamic correla- 
tion length £(tq,). This is qualitatively consistent with 
the fact that the analysis of facilitated models suggests 
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that the length I* is actually different from a dynamic 
heterogeneity length £. For example, for the so-called 
East model, one finds the following relation between the 
dynamic heterogeneity length and the a relaxation time, 
£ ~ exp{-y/m(T Q )} 48J. This relation, for our range of 
correlation lengths and relaxation times, provides as good 
fit to our data as the logarithmic relation £(t q ) ~ ln(r a ) 
discussed above. 



VI. DISCUSSION 

We presented a new computational method to calcu- 
late a dynamic correlation length £(t) characterizing the 
spatial extent of heterogeneous dynamics on time scale t 
and used this new method to obtain a number of results 
pertaining to dynamic heterogeneity. 

Our method combines direct simulational evaluation 
of the four-point structure factor Si(q;t) for non-zero 
wave-vectors with an independent calculation of the dy- 
namic susceptibility Xi(t) that accounts for fluctuations 
suppressed in the simulational ensemble via procedures 
derived by Lebowitz et al. [25] ■ Using the independently 
obtained dynamic susceptibility as the q — > limit of 
Si(q;t) facilitates analyzing the small q behavior of the 
four-point structure factor. We found that an Ornstein- 
Zernicke fits worked well if we restricted our fits such that 
q < 1.5/£(t). This procedure allows one to evaluate the 
dynamic correlation length from simulations of moder- 
ately large systems. We also found that the calculation 
of £(t) from the direct space four-point correlation func- 
tion Gn(r; t) is difficult due to the difficult to account for 
finite size and ensemble dependencies. 

We studied the volume fraction and the time depen- 
dence of the dynamic correlation length. We also ex- 
plored relationships between the length, the dynamic sus- 
ceptibility, and quantities characterizing the average dy- 
namics, the a relaxation time and the self-diffusion coef- 
ficient. 

First, we found that £(r Q ) ~ (0o~ <P)~ 2 provides a good 
description of the data. We note, however, f (T a ) can also 
be fitted by other functions. We also found that mode- 
coupling like power law £(tq,) ~ (4> c — </>)~ 75 provides a 
good description of £(T a ) for the mode-coupling theory 
range of volume fractions, but the exponent 7^ = 0.5 dif- 
fers from the inhomogeneous mode-coupling theory pre- 
diction of 7£ = 0.25 [2D E3- 

Next, we studied the time dependence of Xi{t) an d 
While we did find a power law dependence on time, 
Xi(t) ~ t C i f° r times around the a relaxation time, the 
exponent c was volume fraction dependent and decreased 
with increasing ef>. Surprisingly, we found that for a range 
of times between the /3 and a relaxation times the dy- 
namic correlation length followed a master curve inde- 
pendent of <j) until it reached a volume fraction dependent 
plateau value. The dependence of £(i) on time could be 
fitted with a simple £(t) ~ ln(t) relation. The plateau 
value, £ max , was reached at a characteristic time, r max . 



We found that r max exceeds and grows faster with in- 
creasing volume fraction than the a relaxation time. 

We examined the correlations between £(r a ), r Q , and 
D, and we found that mode-coupling like power law fits 
provide a good description of the data for 0.55 < <f> < 
0.58. We found deviations from these fits as 4> c — 0.59 
is approached. While the mode-coupling exponents for 
T a and D agree reasonably well with the mode-coupling 
predictions, the exponents for £(r Q ) and xa{ t o) do not. 
We also find that x±{ T a) ~ £,{ T a) 3 , which docs not agree 
with the inhomogeneous mode coupling prediction. 

Finally, we found that an exponential dependencies of 
the a relaxation time and the self-diffusion coefficient on 
the dynamic correlation length evaluated at the a relax- 
ation time, T a ~ exp(£(r Q )) and D ~ exp(— £(t q ) ' 6 ), 
describe our data well. This is consistent with the spirit 
of Adam-Gibbs and Random-First-Order- Transition the- 
ories. The values of the scaling exponents are incon- 
sistent with the traditional Adam-Gibbs picture where 
the relaxation, either r Q or l/D, behaves as exp(^). 
Note, however, that Adam-Gibbs and Random-First- 
Order- Transition theories are formulated in terms of the 
characteristic size of dynamically correlated regions, £ s 
whereas we calculated and examined the dynamic cor- 
relation length. Further work is required to clarify the 
connection between these lengths. 
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Appendix A: Characteristic Volume Fractions 

To test theories of the glass transition it is common to 
fit experimental and simulation data to different func- 
tions. The quality of the fits can vary depending on 
the range used for the fits and the proximity to any 
singularity implied by the fitting function. In this sec- 
tion we examine some commonly used fitting functions 
to various <p dependent quantities, a mode-coupling like 
fit (<fc c — 0)~ 7 , a Vogel-Fucher-Tamman (VFT) like form 
exp[A(tj) VF — and a form suggested by Berthier 

and Wittcn (BW), exp[ J B(^ - <^>)" 2 ] G2- The mode- 
coupling fits are used to determine a range of volume 
fractions where the mode-coupling like power laws pro- 
vide a good description of the data. We use this range of 
4> to compare our simulation results to the predictions of 
the mode-coupling theory and the inhomogeneous mode- 
coupling theory. Outside of this range we do not expect 
the mode-coupling theory to provide a very good descrip- 
tion of the dynamics. The goal of this appendix is to 
examine our results and previous arguments in the liter- 
ature to find the best unified description of the data. To 
achieve this goal, we not only examine our data, but also 
use results from earlier investigations {521 149) . 
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TABLE I: Fits to a mode-coupling like power law over differ- 
ent ranges of volume functions. The number in parenthesis 
represents the uncertainty in the last digit. The bottom set 
of data is with <f) c = 0.59 fixed. 



variable 


4>c 


7 


fit range 


D 

{ \ 
Xi\T a ) 


0.5874(4) 
0.5911(7) 
0.58(4) 
0.586(4) 


2.21(6) 
2.13(6) 
1.3(2) 
0.47(9) 


0.55 < (j) < 0.575 
0.55 <(/>< 0.575 
0.55 < <p S 0.575 
0.55 < 4> < 0.575 


T a 

D 

Xi{T a ) 


0.5901(6) 
0.5950(4) 
0.593(3) 
0.61(2) 


2.47(6) 
2.42(5) 
1.7(2) 
0.9(4) 


0.55 < <f) < 0.58 
0.55 < (f) < 0.58 
0.55 <<j)< 0.58 
0.55 < <t> < 0.58 


T a 

D 

Xi{T a ) 


0.59 
0.59 
0.59 
0.59 


2.43(1) 
1.94(3) 
1.46(4) 
0.50(3) 


0.55 < <t> < 0.58 
0.55 < (f) < 0.58 
0.55 <<j)< 0.58 
0.55 < <t> < 0.58 



First we examine mode-coupling like fits to the a re- 
laxation time r a and the self-diffusion coefficient D. The 
mode-coupling theory predicts a power law divergence of 
T a and a power law vanishing of D with the same expo- 
nent 7. It has been found in several numerical studies 
of the mode-coupling theory that 7 « 2.46 [251 [3H]- O ne 
of the difficulties in performing these fits is that, while 
the mode-coupling theory predicts that Stokes-Einstein 
relation is obeyed, r a D — const., |36J, this relation is vi- 
olated in simulations and experiments. Importantly, in 
most simulations the violation of the Stokes-Einstein re- 
lation is apparent already in the regime in which power 
law fits are applicable. Thus, the best one can do is to 
use different exponents for r a and D and the same mode- 
coupling transition volume fraction or temperature. A 
seemingly worse alternative is to force the same expo- 
nent and obtain two different mode coupling transition 
points. 

We fit r Q , 1/D, £(r a ), and Xi{ T a) to power laws 
of the form a((f> c — </>) -7 for 0.55 < <fi < 0.575 and 
0.55 < 4> < 0.58. The results are summarized in Table [I] 
We find that cj) c varies from 0.58 to 0.61, but a consistent 
value is around (f> c « 0.59. Since </> c = 0.59 is consistent 
with our results, has been used in the literature previ- 
ously [27], and coincides with the onset of "hopping" like 
motion observed in Sec. |III| we fix <\> c — 0.59 in this work. 
We also identify 0.55 < <fi < 0.58 as the mode-coupling 
regime, but this should be considered as only an approx- 
imate regime where the mode-coupling theory provides a 
reasonable description of the data. 

Having chosen a common value for the mode-coupling 
transition volume fraction we redo the fits for 1/D, r Q , 
Xi( T a)i and £,( T a) keeping <j) c = 0.59 fixed. The fit pa- 
rameters are the bottom set in Table U 

We now look at the fate of the system beyond the 
mode-coupling regime. To this end we examine the re- 
sults of fits to a VFT like functions and a BW like func- 
tion for r a . We do not fit other variables since we do are 
not sure whether the same functions can be used. Note, 




FIG. 15: Various fits to T a described in the text. The dot- 
ted line is a combination of the fit £(r a ) ~ (4>o — 4>)~ 2 
and the correlation T a ~ exp[fc£(r Q )]. The dashed line is 
the Vogel-Fucher-Tamman fit and the solid line is a fit to 
T a ~ exp[B(0 o - 4>)~ 2 ]- 



however, that £(T a ) is closely tied to r Q , thus we expect 
that fits to £(t q ) result in the same conclusions. 

A fit to ln(T Q ) = ln(7y) + A{<j>v — gives 4>v — 

0.6122±0.0005, t v = 140±7 and A = 0.222±0.005 where 
we fit 0.55 < < 0.5905. Next we fit ln(r Q ) = ln(r ) + 
B(cj) Q -(j))- 2 , which gives t = 456±37, B = 0.017±0.001, 
and 0o = 0.635 ± 0.002. Shown in Fig. [15] are these 
fits along with the results obtained by fitting £,( T a) — 
£0 + C(0o — <f>) an d then using r Q ~ exp(/c£), which we 
refer to as the correlation length fit. The VFT fit is the 
best fit over the largest range of </>, thus one would choose 
this fit based on the fit quality alone. However, the VFT 
fit results in a critical volume fraction that appears to be 
too small when compared with earlier results from the 
literature. 

One result is the dynamic scaling argument of Berthier 
and Witten [32] who found <j) = 0.635 ± 0.005 and 
r a ~ exp[B((j) - (j))- 5 } with 5 = 2.2 ± 0.2. This is in 
remarkable agreement with our fits to r Q and ^(t q ). An- 
other is the work of Odriozola and Berthier [3S] who 
found no evidence of a thermodynamic transition for 
(f> < 0.63 by utilizing a replica exchange Monte Carlo 
algorithm to examine the equation of state. Finally, we 
find that t„w5x 10 8 for a 1 000 particle simulation at 
(j) = 0.6. This value agrees well with the BW fit that 
gives 4.85 x 10 8 for 4> = 0.6, but is orders of magnitude 
different than the prediction of the VFT fit, 1.1 x 10 10 for 
4> = 0.6. While we expect a 1000 particle system to be 
too small for <f> = 0.6 and, in particular, we expect that at 
this volume fraction £(T a ) is larger than half of the 1 000 
particle system simulation cell, this result does provide 
some evidence for the BW fit. However, large, fully equi- 
librated simulations at (j) larger than those utilized in our 
simulations are needed to test the proper functional form 
of the divergence. 
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We should emphasize here that our simulations cover 
both the mode-coupling-like regime and a new regime in 
which the data are consistent with ln(r Q ) ~ (</>o — 4>)~ 2 - 



Appendix B: Calculation of X4(£) and £(t) 

Here we describe the details of the calculation of the 
dynamic susceptibility Xi{t) and the dynamic correlation 
length £(t). The outline of this appendix is as follows. 
First, we present a derivation of Eq. ( 12 ). Next, we show 



that Xiip) calculated from this equation agrees with an 
independent extrapolation of S±(q;t) to q = 0. Finally, 
we use X4 W calculated from formula ( 12 1 as the point at 
Si(q = 0; t) in a fitting procedure to find the most reliable 
result for the dynamic susceptibility and the dynamic 
correlation length. 

In our simulation, the volume fraction and the con- 
centration of particles is fixed. Thus, fluctuations of the 
volume fraction and concentration do not contribute to 
the direct calculation of the dynamic susceptibility. To 
account for these fluctuations we follow the procedure in- 
troduced by Lebowitz et al. |Hj . We start by considering 
an ensemble where the number of particles can fluctuate 
and calculate the first order corrections to Xi{t) calcu- 
lated in an ensemble where the volume fraction and the 
concentration are held constant. Consider 



X*(t) 



(/*1.P2,V) 



5 fewn(t) 5 (l>m(*) 
\ \ n / \ m / 



(Bl) 



where () x denotes an ensemble where X is held constant. 
In Eq. |Bl| the chemical potentials of both small and large 
particles, fj,i and pt 2 , and the volume (V) are held con- 
stant, but the numbers of particles, N± and N 2 , are al- 
lowed to fluctuate. Since the volume is held constant 
in all the ensembles considered, we will not explicitly 
indicate the constant V in what follows. According to 
Eq. 2.11 of Lebowitz et al. gS], 




+2<5iW 2 > Mii/J2 

It Wn /N 1 ,N 2 n 



dN 2 



(B2) 



We utilize the relationship 



lim \/x n x m S nm (q) 

\J x n x m 5 nm , (B3) 



where S nm (q) is the partial structure factor and x n — 
N n /N. We also recognize that En w ™W)jv = (Ni + 
N2)F (t). Finally, we replace the differentiation with re- 
spect to the numbers of particles with the differentiation 
with respect to the volume fraction and the concentration 
and in this way we obtain 



X4.(t) = Xi{t)N + x\ H i + x<t>XcH 2 + X 2 C H 3 

+F (t) 2 H 4 + F (t) X ^H 5 + F (t) Xc H 6 



(B4) 



where Xx = dF (t)/dx. The H n are functions of S nm , 
and are given by 



Hi = 

Hi = 

H 3 = 

H A = 

H 5 = 



( y) [4xiSu + 2d\dl^i- 2 S l2 + d 6 2 x 2 S 22 ] 

(B5) 

~ [d 3 1 x 1 x 2 S 11 - d\x ly /x 1 x 2 Si2 

+d\x 2 ^x\x 2 Si 2 - d\x\x 2 S 22 \ (B6) 
x\x\S\i_ — 2x\x 2 ^J x\x 2 S\ 2 4- x\x 2 S 22 (B7) 
xiSix + 2^/x~[x 2 ~Si 2 + x 2 S 22 (B8) 

^- [dfxiSn + {d\ + dl)^/xix 2 S 12 + d 2 x 2 S 22 \ 

(B9) 



H 6 = 2[xiX 2 Sn + (x 2 - Xi)y/xix 2 Si 2 - XiX 2 S 22 ] . 

(BIO) 

To calculate H n we fit the wave vector dependent ver- 
sion of H ni i.e. expressions (B5|-(B10| with S nm re- 
placed by S nm (q), to a wave- vector independent con- 
stants for q < 0.6. Due to noise in our data we cannot 
perform a more accurate extrapolation. We checked this 
approach by using the same procedure to calculate the 
pressure using the partial structure factors. We checked 
that the pressure obtained from the q — > limit of the 
structure factors agrees with the pressure obtained from 
the extrapolation of the pair correlation function to con- 
tact. 



To verify Eq. (B2) and to check its accuracy, we ex- 
trapolated 5.4 (g; t) to q — > by fitting Si{q;t) obtained 
for non-zero wave-vectors to an Ornstein-Zernicke func- 
tion. We compared the resulting lim 9 _i.o S^{q;t) to Xi{t) 
obtained from Eq. (B2). The extrapolation agreed to 



within error and thus we concluded that Eq. ( B2 ) pro- 



vided a good means to calculate Xi{t)- Subsequently, we 
used Eq. (B2) as the q — value of Si{q;t) in fitting 



procedures. 

It is important to recognize that the above described 
verification of Eq. ( B2 ) requires a rather large system 



size. In particular, we could only perform it using N 
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FIG. 16: The dynamic correlation length £(r a ) obtained using 
the different fitting procedures described in the text. The 
number correspond to the different fitting procedures. 




FIG. 17: Scaling plot S 4 (q; r a ) / S 4 {q 
for the 80 000 particle simulations. 



0;r a ) versus q£{r a ) 



80 000 particles system. Once using Xi{t) obtained from 
Eq. ( |B2[ ) as the q — value of Si{q;t) is accepted, we 
were able to use moderately large systems (N = 10 000 
particles) . 



We fit Si(q;t), using Eq. ((B2J as the q 
S4(q',t), to several functions of the form 



value of 



S 4 (q;t) 



A 

W+W + (i + 



B 



(Bll) 



where all the fitting parameters are time dependent. We 
performed the following fits: (1) set C = and B = 0, i.e. 
an Ornstein-Zernicke type fit; (2) set B=0 which gives a 
function suggested by the inhomogeneous mode-coupling 
theory [3T]; (3) set A = X4(*)U,c and C = 0, which results 
in a function suggested by field theoretic considerations 
[221123]. We also fit S^q; t) to a function utilized by Stein 
and Andersen [13], \u[S 4 (q;t)} = ln(A) - [£q] 2 + Cq 4 , 
procedure (4). All of the fits results except for proce- 



dure (3) results in statistically the same length, Fig. 16 
if we restrict the fits as follows. For procedure (1), the 
Ornstein-Zernicke fits, we only fit to q < 1.5/£ and for 
the fit to \n[S4(q;t)] we only fit q < l/£. Procedure (3) 
resulted in an £ approximately 1.2 times smaller than the 
other procedures at every volume fraction, thus none of 
the conclusions of this work changes due to utilizing that 
fitting function. For volume fractions beyond our abil- 
ity to study, it may be found that £ determined through 
procedures (1), (2), and (4) is not simply a factor of £ 
found using procedure (3). As a final check, we used £ ob- 
tained from the Ornstein-Zernicke fit to check the qual- 
ity of overlap of S±(q;Ta)/ S±(q — > 0;r Q ) versus q£(T a ), 
Fig. [TTJ and we find the overlap to be very good. The 
results shown in Figs, ph, (111]), (p), ph, (111), and 



(17) are found by the Ornstein-Zernicke fits 
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